source("/users/brucedesmarais/Dropbox/SenPE/Code/CreateSenPENets.R")

library(igraph)
community <- list()

for(t in 1:length(mats)){
mati <- mats[[t]]
ntype <- c(rep(0,nrow(mati)),rep(1,ncol(mati)))
el <- NULL
for(i in 1:nrow(mati)){
	for(j in 1:ncol(mati)){
		if(mati[i,j]==1){
			n1 <- i-1
			n2 <- nrow(mati)+j-1
			el <- rbind(el,c(n1,n2))
			}
		}
	}

# Project the bipartite network

g <- graph.bipartite(ntype,c(t(el)))
g2 <- bipartite.projection(g)[[2]]
m <- as.matrix(g2)
edge_g2 <- cbind(m[4][[1]],m[3][[1]])
wg2 <- numeric(nrow(edge_g2))
for(i in 1:nrow(edge_g2)){
		indi <- edge_g2[i,1]
		indj <- edge_g2[i,2]
		wg2[i] <- sum(mati[,indi+1]*mati[,indj+1])
	} 

cg2 <- fastgreedy.community(g2)
cg2w <- fastgreedy.community(g2,weights=wg2)
mem <- community.to.membership(g2,cg2$merges,which.max(cg2$modularity))
memw <- community.to.membership(g2,cg2w$merges,which.max(cg2w$modularity))
nett <- get.adjacency(g2)
nett[edge_g2+1] <- wg2
community[[t]] <- list(cg2,cg2w,mem,memw,nett)
}

library(foreign)
senData <- read.dta("/Users/brucedesmarais/Dropbox/SenPE/Data/fowlersenpenbinomdata.dta")

procName <- function(x){
	xvec <- character(length(x))
	for(i in 1:length(x)){
		chari <- x[i]
		comind <- which(strsplit(chari,split="")[[1]]==",")
		if(length(comind)==1){
		if(comind +1 < nchar(chari)){
			chari <- substr(chari,1,comind+2)
			}
			}
		xvec[i] <- chari
		}
		xvec
	}
	
write.dta(senData)

fow_cong <- paste(senData$labels,senData$congress,sep="")
pe_cong <- paste(procName(id_master$sen_pe),id_master$Congress,sep="")
fow_master <- paste(id_master$cosp_summary,id_master$Congress,sep="")
sen_cong <- pe_cong[match(fow_cong,fow_master)]
comm <- rep(10,length(sen_cong))
for(i in 5:13){
	comi <- community[[i]][[4]]$membership
	memi <- procName(colnames(mats[[i]]))
	cong <- congs[[i]]
	mem_cong <- paste(memi,cong,sep="")
	cind <- match(mem_cong,sen_cong)
	comm[cind[!is.na(cind)]] <- comi[!is.na(cind)]
	print(mean(is.na(cind)))
	}
	
#senData <- data.frame(senData,comm)
write.dta(senData,"senDataCommunity.dta")

distData <- NULL
for(c in 97:105){
	subDat <- subset(senData,senData$congress==c)
	distDatc <- NULL
for(i in 1:ncol(subDat)){
	distDatc <- cbind(distDatc,c(dist(subDat[,i])))
	}
	distData <- rbind(distData,cbind(distDatc,c))
	print(c)
	}

smCmte <- NULL
smComm <- NULL
smCong <- NULL
smPrty <- NULL
dsIdeo <- NULL
dsSenr <- NULL
cmtes <- paste("cmte",1:23,sep="")
for(c in 97:105){
	subDat <- subset(senData,senData$congress==c)
for(i in 2:nrow(subDat)){
	for(j in 1:(i-1)){
	cmte1 <- t(round(subDat[i,cmtes],-2))
	cmte2 <- t(round(subDat[j,cmtes],-2))
	
	nCmte <- sum(is.element(na.omit(unique(cmte1)),cmte2))
	sComm <- subDat$comm[i] == subDat$comm[j]
	if(sum(c(subDat$comm[i],subDat$comm[j])==10)>0) sComm <- NA
	smPty <- subDat$party[i]==subDat$party[j]
	smCmte <- c(smCmte,nCmte)
	smComm <- c(smComm,sComm)
	smCong <- c(smCong,c)
	smPrty <- c(smPrty,smPty)
	dsIdeo <- c(dsIdeo,abs(subDat$ideol1[i]-subDat$ideol1[j]))
	dsSenr <- c(dsSenr,abs(subDat$senior[i]-subDat$senior[j]))
	}
	}
	print(c)
	}

prTab	 <- matrix(0,9,4)
frTab	 <- matrix(0,9,4)
nTab	 <- matrix(0,9,4)
lcTab	 <- matrix(0,9,4)
ucTab	 <- matrix(0,9,4)
for(c in 97:105){
	for(sm in 0:3){
		prTab[c-96,sm+1] <- mean(subset(smComm,smCong==c & smCmte==sm),na.rm=T)
		frTab[c-96,sm+1] <- sum(subset(smComm,smCong==c & smCmte==sm),na.rm=T)
		nTab[c-96,sm+1] <- length(na.omit(subset(smComm,smCong==c & smCmte==sm)))
		bt <- binom.test(sum(subset(smComm,smCong==c & smCmte==sm),na.rm=T),length(na.omit(subset(smComm,smCong==c & smCmte==sm))))
		lcTab[c-96,sm+1] <- bt$conf.int[1]
		ucTab[c-96,sm+1] <- bt$conf.int[2]
		}
	}

library(gplots)

lciCmte <- numeric(9)
lciPrty <- numeric(9)
uciCmte <- numeric(9)
uciPrty <- numeric(9)
dprCmte <- numeric(9)
dprPrty <- numeric(9)
for(c in 97:105){
		smCmteC <- as.numeric(subset(smCmte,smCong==c)>0)
		smCommC <- subset(smComm,smCong==c)
		smPrtyC <- subset(smPrty,smCong==c)
		tabPrt <- table(smPrtyC,smCommC,useNA="no")
		tabCmt <- table(smCmteC,smCommC,useNA="no")
		prtPrt <- prop.test(tabPrt)
		prtCmt <- prop.test(tabCmt)
		dprCmte[c-96] <- -diff(prtCmt$estimate)
		dprPrty[c-96] <- -diff(prtPrt$estimate)
		lciCmte[c-96] <- prtCmt$conf.int[1]
		uciCmte[c-96] <- prtCmt$conf.int[2]
		lciPrty[c-96] <- prtPrt$conf.int[1]
		uciPrty[c-96] <- prtPrt$conf.int[2]
	}


par(mfrow=c(3,3))
i <- 1

pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/commCmte.pdf",height=3.5, width=6, family="Times", pointsize=15)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plotCI(x=97:105,y= dprCmte,li=lciCmte,ui=uciCmte,col="blue",xaxt="n",ylab="Difference in Proportion",xlab="")
text(97:105, par("usr")[3] - 0.015, srt = 45, adj = 1,
    labels=97:105, xpd = TRUE,lwd=2.5)
abline(h=0,lty=2,lwd=1.5)
axis(1,labels=F,at=97:105)
dev.off()

pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/commPrty.pdf",height=3.5, width=6, family="Times", pointsize=15)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plotCI(x=97:105,y= dprPrty,li=lciPrty,ui=uciPrty,col="blue",xaxt="n",ylab="Difference in Proportion",xlab="",gap=.5)
text(97:105, par("usr")[3] - 0.05, srt = 45, adj = 1,
    labels=97:105, xpd = TRUE,lwd=2.5)
abline(h=0,lty=2,lwd=1.5)
axis(1,labels=F,at=97:105)
dev.off()

pcpPrty <- NULL
pcpCmte <- NULL
pcpBoth <- NULL
for(c in 97:105){
	both <- glm(smComm~smPrty+smCmte,subset=smCong == c & ! is.na(smCmte) & !is.na(smPrty),family=binomial)
	prty <- glm(smComm~smPrty,subset=smCong == c & ! is.na(smCmte) & !is.na(smPrty),family=binomial)
	cmte <- glm(smComm~smCmte,subset=smCong == c & ! is.na(smCmte) & !is.na(smPrty),family=binomial)
	prBoth <- as.numeric(predict(both)>0)
	prCmte <- as.numeric(predict(cmte)>0)
	prPrty <- as.numeric(predict(prty)>0)
	y <- both$y
	pcpPrty <- c(pcpPrty,AIC(prty))
	pcpCmte <- c(pcpCmte,AIC(cmte))
	pcpBoth <- c(pcpBoth,AIC(both))
	}
	 
